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FOREWORD 

This  work  was  completed  on  the  DIDS  project  of  the  U.   S. 
Army  Strategic  Communications  Command,   Fort  Huachuca,  Arizona 
85613  as  part  of  the  work  of  interagency  work  order  DAHC  -20-69-C  -Oil  6. 

The  computer  programs  were  developed  for  the  Defense  Atomic 
Support  Agency  under  interagency  work  order  805-69o     Test  computa- 
tions made  in  connection  with  this  development  are  also  included. 
Permission  to  quote  data  points  plotted  in  figures  30,    31,    32  and  33 
(Morgan,    1966)  was  granted  by  the  Department  of  the  Army,   Office 
of  the  Secretary  of  the  Army,   Office  of  Civil  Defense. 

The  FORTRAN  listing  of  the  computer  program  described  in 
Section  7  is  available  from  the  authors  upon  request. 
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ABSTRACT 

The  propagation  of  LF  and  VLF  terrestrial  radio  waves  is 
simulated  on  the  computer  with  a  precise  theory  of  spherical  wave 
functions  of  integer  order.     Anisotropic  reflections  from  an  electron 
and  ion  density  profile  of  the  ionosphere  are  accommodated  by  the 
analysis.     Variations  in  the  profile  and  the  magnetic  field  'with  distance 
along  the  propagation  path  can  also  be  introduced  ad  hoc. 

The  analysis  method  is  incorporated  into  a  computer 
program  which  is  written  in  a  flexible  manner  so  that  it  can  be 
applied  to  a  variety  of  scientific  studies  of  low-frequency  ionospheric 
wave  propagation.     Calculations  of  current  interest  are  presented  to 
illustrate  use  of  the  technique.     Amplitude  perturbations  of  approximate- 
ly 1  5  decibels  are  described  at  60  kHz  as  a  result  of  an  electron  density 
profile  perturbation.     It  is  also  found  that  the  26  kHz  signal  is  compara- 
tively insensitive  to  a  variation  in  magnetic  dip  angle  when  propagation 
occurs  in  the  daytime  near  the  magnetic  equator.     Thus,   using  values 
of  the  magnetic  field  vector  at  the  center  of  the  propagation  path  is  a 
good  approximation  in  certain  circumstances. 
KeyWords:         MF,    LF,    VLF;    low  frequency  wave  propagation; 

propagation;    terrestrial  radio  waves;    low  frequency 
radio  waves;    radio  waves;    spherical  waves; 
spherical  wave  functions. 
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1.     INTRODUCTION 

The  formulation  of  a  theory  of  propagation  for  MF,    LF,   and 
VLF  ionospheric  radio  waves  using  spherical  wave  functions  of  integer 
order  as  building  blocks  was  presented  in  a  previous  work  (Johler, 
1970).     Ordinarily,    the  spherical  wave  functions  of  integer  order  are 
found  to  be  useful  only  at  ELF  or  for  propagation  problems  in  which    ka 
is  small,   where    a    is  the  radius  of  the  sphere  and    k    is  the  wave 
number  (Johler  and  Lewis,    1969). 

In  the  approach  presented  in  this  paper,    this  restriction  is  no 
longer  necessary  because  of  the  properties  of  the  j -series  summation 
used.     Indeed,   numerical  testing  of  the  method  at  medium  frequencies 
(MF)  indicates  that  the  upper  frequency  limit  of  usefulness  of  spherical 
wave  functions  had  not  as  yet  been  found.     The  theory  has  been  trans- 
lated into  a  FORTRAN  computer  program  for  use  in  scientific  investi- 
gations of  the  ionospheric  radio  waves.     The  philosophy  in  the  computer 
program  design  emphasized  flexibility  in  the  application  of  the  analysis 
method  to  a  variety  of  scientific  studies  concerned  with  MF,    LF  and 
VLF  wave  propagation.     Obviously,    extensive  calculation  of  the  field, 
such  as  that  employed  in  system  studies,   would  require  optimization 
of  the  computer  program  for  a  specific  system  application  to  insure 
economy  of  computer  time. 


This  paper  demonstrates  the  application  of  the  method  of  spheri- 
cal wave  functions  of  integer  order  to  LF  and  VLF.     Calculations  of 
interest  are  presented  and  compared  with  experimental  data.     The 
FORTRAN  computer  programs  necessary  for  the  calculations  pre- 
sented are  also  described. 

2.     COMPUTATION  AND  DISCUSSION 

The  theoretical  approach  used  in  this  paper  is  the  geometric 
series  propagation  formula  (Johler,    1970) 

CO 
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where    Er    0    is  the  ground  wave  (Johler  et  al.  ,    1956;   Johler  and 
Morgenstern,    1965;    Johler,    1969a, b)    and    E      i    are  the  terms  of  the 
geometric  series,     j  =  1,    2,    3  •  *  •    ,    depicted  as  optical  rays  in 
figure  1.     It  should  be  emphasized  that  there  is  no  set  of  optical  rays  of 
the  homogeneous,   plane  wave  type  that  will  give  the  exact  solution  pre- 
sented in  this  analysis  except  as  a  very  approximate  estimate.     How- 
ever,  a  rigorous  set  of  n-waves,    discussed  in  detail  in  section  3  will 
be  referenced  in  the  physical  discussion  to  follow.     The  n-waves  are 
almost  equivalent  to  homogeneous  plane  waves  with  respect  to  the 
ionosphere  boundaries. 

If  the  terrestrial  waveguide  is  assumed  to  be  sharply  bounded 
at  the  top  with  a  reflection  coefficient,     T,  #   =   -1,   various  ionospheric 

waves,     Erj  j     ,     j  =  1,    2,    3  *  *  "    ,    can  be  calculated  (Johler,    1970)  from 

t 
the  formula    , 


tQuantities    B,G,R,  ,p0   etc.  ,   are  defined  in  section  6. 
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Each  term  of  this  summation  (2)  is  a  spherical  wave  of  integer  order  or 
an  n-wave  (Johler,    1970).     The  theoretical  basis  of  such  n-waves  is 
discussed  in  section  3.     The  spherical  reflection  coefficient  for  a 
TM-wave  (vertically  polarized  electric  vector)  on  the  ground,     R,  ,     is 
a  function  of  the  ground  conductivity,     o,    mhos/m,    as  is  the  wave 
number,     k_2  .     If    o  =  5  mhos/m  (sea  water),    the  amplitude,    |Er    i  \, 
is  illustrated  in  figure  Z,   j  =1  through  10.     The  height  of  the  waveguide 
is    h  =  65  km.     The  frequency,    f  =  10  kHz.     The  curves  presented  are 
therefore  independent  of  the  complicated  ionosphere  reflection  coeffi- 
cient and  (except  at  very  short  distances)  illustrate  the  basic  character 
of  the  various  terms  of  the  j -series  (or  hop  series).     Each  term,     Er,j, 
at  short  distance,    d/ j,    is  a  monotonic  curve.     Undulations  then  appear 
as  a  function  of  distance  as    d/  j    is  increased.     In  the  region  of  the 
geometric -optical  horizon  of  the  ray  theory,   and  beyond,   as  depicted 
in  figure  1,   the  wave  theory  gives  a  curve  that  shows  an  exponential 
type  of  attenuation  characteristic  of  diffracted  waves.     The  standing 
wave  as  a  function  of  distance  in  the  intermediate  region  is  evident  at 
a  certain  critical  angle  of  incidence,     cpi  ,     on  the  ground.     This  sug- 
gests lateral  wave  displacement  (Brekhovskikh,    I960),    since  a  critical 
domain  of  the  complex  angle  of  incidence  on  the  ground  or  the  iono- 
sphere is   necessary  to  excite  such  waves  (Johler,    1970).     Since, 
T«,    =   -1»   these  are  lateral  waves  along  the  ground. 

At  great  distance,    the  waves  diffract  around  the  curvature  of 
the  ground.     In  such  a  region  a  considerable  number  of  n-waves  or 
terms  of  the  summation  over  n,    discussed  in  section  3,   are  important 
in  the  series  (2).     As  a  consequence  of  this  fact,    the  curves  are  not 
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precisely  exponential  (linear  on  a  semilog  graph).     In  figure  3,    the 
ground  wave,     Er   0  ,     is  also  shown  for  comparison.     A  higher  frequency 
(2  6  kHz)  is  illustrated  in  figure  4,   assuming  a  ground  conductivity 
0  =  0.  005  mhos/  m. 

At  short  distances  the  j -series  terms  ordinarily  have  small 
amplitude.     With    T,.   =   -1,    the  amplitudes  illustrated  in  figures  2,    3 
and  4  are  abnormally  large.     Thus,    except  at  very  short  distances,   the 
basic  features  shown  in  the  figures  2,    3,    4  dominate  the  terms  of  the 
j -series.     More  complicated  effects  can  be  introduced  by  using  more 
realistic  reflection  coefficients  for  the  ionosphere.     An  exponential 
model  of  the  ionosphere  and  a  two  layer  model  ionosphere  developed  by 
Thrane*  in  1969  are  illustrated  in  figure  5„     An  alternatie  "two -layer 
model"  profile  developed  by  Bain  and  May  (1967)  is  illustrated  in 
figure   6.     In  a  "two -layer  model"  the  electron  density  tends  to  form  a 
clear  maximum  in  the  D-region  near  60  km.     The  electron  density  then 
decreases  before  increasing  again.     Both  electron  density  and    d-c 
conductivity  for  the  electrons,    o~#    , 

Ne2  ,~a 

a     =  (3) 

e        meV 

are  shown  in  the  figures  as  a  function  of  altitude  above  the  ground  level, 
km.     A  theoretical  daytime  noon  model  developed  by  W.  S.  Knapp**  in 
1967  based  on  computer  programs  of  the  type  given  by  D.  C.  Hyovalti*** 
in  1965  ,     is  shown  in  figure  7.     Both  the  electron  density  and  the 
positive  ion  density  are  given.     The  negative  ion  density  can  be  obtained 
from  conservation  of  change. 

N+    =  Ne  +  N-     .  (4) 


*      E.  V.    Thrane,   Private  Communication,   ITS,   ESSA,    Boulder,    Colo. 
*#   W.  S.   Knapp,   Private  Communication,   G.E.  Tempo,   Santa  Barbara, 

Calif. 
#**D.  C.  Hyovalti,  Unpublished  Report,    ITS,   ESSA,    Boulder,   Colo. 
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Then,   the  total    d-c  conductivity  is  given  by 

0.   ,  =  H^L    +  2  2k*  ,  (5) 

where 

v'       =  v         Els (a  =  e,+) 

a'°  a'°    ma  +  m0     '    v 

and  where  the  mass  m,   of  oxygen  and  the  mass  m.    of  the  electron  can  be 
used.     Of  course,   the  effect  of  ions  is  quite  small  in  the  normal  ionosphere 
since 

me  <  <  m+      .  (6) 

Collision  frequencies  used  by  Johler  and  Berry  (1965)  were  used  in  all 
calculations  involving    v#J  0  .     Positive  ion  collision  frequencies,    v+    0    , 
and  negative  ion  collision  frequencies  are  also  given  by  Hirsch  (1967). 
Under  disturbed  conditions  such  as  polar  cap  events  or  nuclear  events 
(Johler  and  Berry,    1965),    it  may  be  necessary  to  consider  conductivity 
caused  by  ions. 

A  daytime  profile  similar  to  the  profile  in  figure  7  is  given  in 
figure  8.     However,    the  electron  densities  at  the  lower  altitudes  between 
40  and  70  km  become  negligible  more  abruptly  with  decreasing  altitude. 
There  is  some  uncertainty  concerning  the  precise  form  of  the  lower 
boundary.     This  can  affect  the  propagation  of  terrestrial  radio  waves 
in  the  waveguide  below  the  ionosphere.     Belrose  and  Segal  (1967)  give 
comparative  profiles  of  electron  number  density  at  the  lower  boundary 
of  the  ionosphere. 

The  model  shown  in  figure  8  is  compared  in  figure  9  with  a 
nighttime  model  developed  from  similar  computer  calculations.     Also, 
the  corresponding  collision  frequency  is  given.     Of  course,    other  pro- 
files are  available  in  the  literature  that  can  be  used  to  model  the 
ionosphere  (Belrose,    1968).     To  demonstrate  the  use  of  spherical  wave 
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function  methods  we  shall  use  the  models  given  in  figures   5  through  9 
in  the  computer  program  discussed  in  section  6. 

The  models  given  in  figures  7  and  8  are  regarded  by  us  as  more 
representative  of  daytime  noon  conditions  of  the  ionosphere.     The  lat- 
ter is  a  more  sharply  bounded  model.     The  models  given  in  figures  5 
and  6  are  special  models  that  may  be  used  to  explain  special  pheno- 
mena,   such  as  the  effect  of  a  tendency  to  form  a  lower  layer  of 
electrons.     In  such  discussions  the  exponential  model  given  in  figure  5 
is  often  referenced.     The  available  data  on  the  nighttime  ionosphere  is 
not  as  abundant  as  the  daytime  noon  model.     The  nighttime  model 
given  in  figure  9  may  or  may  not  be  typical. 

In  (2),   let 

p,  =  (-1)*  pi  Rl"1        ;  (?) 

then,   an  isotropic  reflection  coefficient  can  be  introduced  by  replacing 
-1     with    Te-J 

Cj   =  (TeJ*  ^0  Rf1  ,  (8) 

where    T09     is  the  reflection  coefficient  for  a  wave  incident  on  the 
ionosphere  with  the  electric  vector  of  the  incident  and  reflected  wave  in 
the  plane  of  incidence.     This  is  called  a  TM-wave  (transverse  magnetic) 
reflection.     If  the  ionosphere  is  anisotropic,     C3     is  more  complicated. 
Thus,    (Johler,    1961), 


Ci  -  T9e 

>g      =    XVe       J.  ee  i      £\.m       X  6m       J.me 


C2   -  Re  Tee    +  Rra   T-m   T, 


2      rp        3 

ee 


C3   -  2Re  Rn  T9e  Tem  Tme  +  R0     T, 

(9) 
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A  simple  system  for  working  out  the  first  few    Cj     is  depicted 
in  figure  10  using  a  geometric  notion  of  rays  reflected  to  and  from 
between  the  two  waveguide  boundaries.     The  rigorous  matrix  form  of 
this  solution  is  given  in  other  papers,   Johler  (1970)  and  Johler  (1966). 
Recently,    Berry,    et  ah    (1969)  have  worked  out  this  matrix  expansion 
so  that  ionosphere  reflection  coefficient  effects  can  be  separated  from 
ground  reflection  coefficient  effects  as  was  originally  postulated  in 
previous  works  (Berry,    1964;   Berry  and  Chrisman,    1965).     It  is 
quite  apparent  from  figure  10  that  the  expressions  for  the    Cj  '  s 
can  become  quite  complicated  and  uninteresting.     However,   the  first 
few  terms  of  the  geometric  series  are  instructive  from  a  physical  point 
of  view.     Thus,    considering  the  first  term,     j  =  1,     of  the  geometric 
series,   the    incident  vertical  electric  polarization  gives  rise  only  to 
TM -wave  propagation  between  the  transmitter  and  the  receiver.     While 
it  is  true  that  an  abnormal  component  is  generated  at  the  ionosphere  by 
the  conversion  coefficient,     Tem,     the  assumed  vertical  polarization  of 
the  source  and  observer  dipoles  does  not  involve  this  component  in  the 
waveguide  propagation  mechanism  shown.     For  two  reflections,     j  =  2, 
the  reflected  TE -wave  (dashed  ray)  is  reflected  at  the  ground  and  then 
converted  back  into  a  TM-wave  at  the  second  reflection  by  the  operation 
Tem  Tne.     Of  course,   the  conversion  coefficient,  Tem ,  indicates  conver- 
sion occurs  in  a  region    [2,  1]  on  the  second  term    j  =  2,   first  reflection 
k  =  1     of  the  set    [j,k],    while    Tn0     indicates  conversion  occurs  in  a 
region  \_2,  2]    at  a  greater  distance  from  the  transmitter.     Obviously  the 
process  can  be  continued  ad  infinitum,   as  indicated  by    j  =  3  and  4  in 
figure  10.     Thus,   in  general  a  region  labeled    [j,k]     is  the    kth   reflec- 
tion (beginning  at  the  transmitter  and  numbering  the  regions  with  the 
k  =  1 ,    2,    3  •  •  ■   integers)  of    Er>  j  .       When    j  =  3    or  greater,   the  TE  - 
wave  reflection  coefficient,     Tn  n  ,   becomes  significant  to  propagation  in 
which  both  transmitter  and  receiver  are  vertical  dipoles.     This  is  a 


consequence  of  the  conversion  effect  of  the  terrestrial  magnetic  field  at 
the  ionosphere  since  the  source  does  not  give  rise  to    TE  -waves.     Thus, 
TM    and    TE    propagation  modes  become  coupled  at  the  ionosphere.    As 
the  terrestrial  magnetic  field  is  theoretically  forced  to  vanish, 
Ten  =  TB  e   =  0,  and  propagation  is  controlled  by    Te,   (TM-waves)  exclu- 
sively.    Note  also  that  the  ground  reflection  coefficient,     Rm  ,     for    TE  - 
waves  is  not  used  in  the  propagation  mechanism  if  the  ionosphere  is 
isotropic. 

It  is  clear  from  the  indexing    [j,k]     that  regions  of  the  ionosphere 
have  been  identified  along  the  propagation  path  between  transmitter  and 
receiver.     These  regions  can  be  tagged  with  the  aid  of  points  represent- 
ing the  center  of  the  reflecting  region    [j,k]    by  drawing  geometric - 
optical  rays  depicted  in  figure  1.     The  waveguide  height  has  been  postu- 
lated as    h.     This  is  a  physical  parameter.     In  actuality,   the  geometric 
series  expansion  is  ordinarily  performed  in  a  thick  bottom  layer  repre- 
senting the  earth-ionosphere  waveguide  as  is  illustrated  in  figure  11. 
In  a  rigorous  solution,    the  ionosphere  is  then  represented  by  concentric 
spherical  shells  to  great  height  such  that  decreasing  the  shell  thickness 
or  increasing  the  number  of  such  shells  or  decreasing  the  height  to  which 
such  shells  are  emplaced  produces  negligible  effect  on  the  field  propagated 
from    S    to    O.     The  precise  thickness  of  the  waveguide  shell  in  which  the 
j -series  expansion  is  performed  is  then  primarily  a  matter  of  phase 
reference  for  the  reflection  coefficient  calculation,   although  a  choice 
can  often  be  found  in  which  the  phase  of  the  reflection  coefficient,     Tc0 
varies  slowly  in  the  summation  process  (2).     This  slowly  varying  phase 
reference  can  be  used  to  improve  the  efficiency  of  the  calculation,    since 
the  number  of  values  of  reflection  coefficient  calculated  can  be  decreased. 
Values  of    Te  ,   =  T,e(n)    are  found  for  each    n    by  interpolation  in  a 
calculated  set  of    Te#(cp1)     as  a  function  of    cpt  ,     related  to    n    by  (see 
sec.    3), 
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where    g  =  a  +  h,     k0   =  — ,  and    a    is  the  radius  of  the  terrestrial  sphere, 
a~  63  60  km.     Although  the  computer  program  can  accommodate  complex 
angles  of  incidence,     cp4  ,     it  is  evident  from  the  approximation  (11)  that 
the  angle  is  almost  real. 

In  the  presence  of  magneto -ionic  coupling  of  the  wave  into  the 
terrestrial  magnetic  field  lines,    the  phase  of  the  reflection  coefficient 
does  not  usually  display  such  a  minimum  phase  variation  with  angle  of 
incidence,     cpt  .     Under  such  conditions  a  quite  high  density  of  values  of 
reflection  coefficient  as  a  function  of    cpt     or    n    must  be  calculated  for 
high  accuracy. 

Ordinarily,   the  bottom  of  the  lowest  spherical  shell  (figs.  1 
through  11)  at  a  radial  distance    r  =  a  +  h,     is  used  as  the  reference 
point  for  the  phase  of  the  reflection  coefficient.     If  a  criterion  of  reflec- 
tion height  is  used  in  such  a  way  that  the  phase  is  slowly  varying  as  a 
function  of    cp±  ,     the  phase  of  the  reflection  coefficient  should  be  re- 
referenced  to  the  height  used  (Budden  1961).     Thus  if    «    means     "is 
replaced  by" , 

T,t  «  T..   exp  [2ik0   coscpi  (h  -  h0)]   ,  (12) 

where  h0  is  the  bottom  of  the  ionosphere  where  the  reflection  coeffi- 
cient is  calculated,  and  h  is  the  height  to  which  the  reflection  coeffi- 
cient is  to  be  re -referenced.  It  is  quite  apparant  that  the  phase  of  the 
reflected  wave  is  retarded  if  h  >  h0  ,  and  advanced  if  h  <  h0  .  Equa- 
tion (12)  is  restricted  in  usage  to  plane  wave  reflection  from  a  plane 


ionosphere.     The  use  of  spherical  shells  would  require  more  compli- 
cated modification  (Johler,    1970).     In  many  calculations  the  use  of  avail- 
able reflection  coefficients  that  are  based  on  plane  wave,   plane  iono- 
sphere models  is  quite  accurate  and  convenient.     A  computer  program 
is  described  in  section  7  for  computing  the  reflection  coefficient  of  an 
anisotropic  ionosphere  using  an  exact  method  for  the  reflection  of  a 
plane -wave  from  a  plane -ionosphere  with  continuous  stratification  of 
the  profile  with  altitude.     It  is  used  with  calculations  performed  with  the 
spherical  wave  function  programs  which  are  also  given  in  section  7. 

Some  further  discussion  of  the  geometric  series  expansion  in 
the  thick  bottom  isotropic  shell,    (fig.    11)  containing  the  transmitter    S 
and  the  receiver    O  (the  source  and  observer)  is  appropriate  to  a  further 
understanding  of  the  nature  of  the  reflection  coefficient.     After  the  height 
of  the  bottom  shell,     h,     has  been  established,    the    jth  term  of  the  series 
(1)  can  be  calculated  by  the  summation  (2).     Each  term  of  the  summa- 
tion (2)  is  a  spherical  wave  of  integer  order  (Johler,    1970)  or  an  n-wave. 
The  theoretical  basis  of  such  n-waves  is  discussed  in  section  3.     One 
such  n-wave  is  depicted  in  figure  11,    for    j  =  1     and    j  =  2.     The  upgoing 
wave  is  shown  with  a  solid  line.     The  downgoing  wave  is  shown  with  a 
large  dashed  line.     Since  a  particular  n-wave  has  been  chosen,    the  angle 
of  incidence  with  the  ionosphere  is  given  by  (10)  or  approximately  by  (11). 
A  short-dashed  line  shows  the  angle  of  incidence  of  the  geometric - 
optical  ray.     Note  that  the  angle  of  incidence  of  the  geometric -optical 
ray  and  the  n-wave  do  not  necessarily  correspond.     One  can  make  the 
conjecture  that  important  contributions  to  the  sum  of  the  n-waves  (2) 
will  occur  when  the  angle  of  incidence  of  the  n-waves  approaches  the 
geometric -optical  wave.     This,    indeed,   has  been  found  to  be  precisely 
the  case  when  the  series  is  summed  for  models  presented  in  this  paper. 
In  fact,    the  n-waves  with  very  steep  angles  of  incidence,     cpt  ,     on  the 
ionosphere  were  found  to  be  severely  attenuated  as  were  the  almost 
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grazing  incidence  waves.     In  the  latter  case,    the  angle  of  incidence 
became  complex    (-Imsincpi    >  >  0) .     Suppose  the  n-wave  depicted  in 
figure  1 1  is  the  lower  limit  n-wave  below  which  smaller    n    values 
given  negligible  contribution.     Then  as  equation  (2)  is  summed,    it  can 
be  imagined  that  the  regions  \.\,  l~\,   [2,1],    [.2,2]    ■  •  •  ,   are  the  important 
regions  for  the  terms  of  the  j  -series.     Then,    instead  of  an  illuminated 
point  used  in  geometric -optical  theory,    we  have  an  illuminated  point  set 
extending  over  a  finite  region  of  the  propagation  path.     It  is  also  clear 
that  the  n-waves  of  small  order    n    must  exhibit  considerable  lateral 
displacement  along  the  ionosphere.     One  must  at  the  same  time  admit 
lateral  displacement  along  the  ground.     This  suggests  physical  modifi- 
cations of  the  analysis  which  can  take  into  account  local  regions  of  the 
ionosphere    [1,1],    [2,1],    [2,2]   *••     and  local  regions  of  the  ground 
[1,1],   [1,2],   [2,1],   [2,2],    [2,3]  •••    ,     with  different  reflection 
coefficients,     T,,(j,k),   R,(j,k)   •••    .     The  reflection  coefficients 
T,,(j,k)    may  be  a  consequence  of  profile  or  magnetic  field  variations 
along  the  propagation  path.     The  reflection  coefficient    R,  (j,k)    may  be 
a  consequence  of  ground  conductivity  or  ground  impedance  changes 
along  the  propagation  path. 

Suppose  for  example,    the  ionosphere  is  disturbed  over  a  region 
of  radius    R,   as  depicted  in  figure  12.     If  the  region    R    were  of  suffi- 
cient size,    it  could  be  represented  by  a  spherical  reflection  coefficient 
T,,     and  a  reference  height  that  may  be  the  reflection  height,     h0    . 
Also,   a  disturbance  at  some  shorter  distance  (fig.    13)  may  have  no 
effect  on  the  term    Er  fl  of  the  geometric  series  and  the  normal  reflec- 
tion coefficient  would  be  appropriate.     But  the    Er  j2    term  of  the  series 
would  require  a  special  reflection  coefficient  and  reflection  height. 
The    Er    3    wave  may  be  unaffected.     But  the  total  field,     Er  ,    of 
equation  (2)  will  be  affected  by  this  disturbance.     The  case  depicted  in 
figure  12  may  also  affect  the    Er  3    term  of  the  geometric  series,   as 


11- 


shown  in  figure  14  to  the  exclusion  of  the    Er    2    term.     Finally,   we  note 
the    Er    4    term  could  also  be  disturbed  by  a  disturbance  in  the  region 
[4,2],    i.e.,    the  fourth  hop,    second  reflection  (fig.    15).     The  flexibility 
described  above  has  been  incorporated  into  the  computer  program 
described  in  section  7„     With  some  imagination,    it  is  not  difficult  to 
envision  further  research  on  the  nature  of  propagation  via  the  D -region 
using  these  computer  programs  as  a  tool.     It  should  be  emphasized  that  a 
sufficient  number  of  n-waves  should  be  included  to  assure  convergence 
of  the  n-series.     Also,   a  sufficient  number  of  j -waves  or  terms  of  the 
j  series  must  be  included  for  the  desired  accuracy.     Usually     3    or 
4  j -waves  will  suffice.     If  a  large  number  are  required,   the  above 
described  additional  flexibility  feature  becomes  difficult  to  apply  for 
obvious  reasons.     Finally,   the  lower  limit  of  the  n-series  summation 
should  be  carefully  determined  with  appropriate  tests  to  be  sure  that 
the  steeply  incident  waves  are  completely  negligible. 

A  rigorous  method  for  calculating  plane -wave,   plane -ionosphere 
reflections  in  the  presence  of  arbitrary  magnetic  induction  and  arbitrary 
variation  of  ion  or  electron  number  density  and  collision  frequency  with 
altitude  was  given  by  Johler  and  Harper  (19  62).     The  use  of  this  method 
was  reviewed  more  recently  by  Johler  (1967).     Thus,    in  addition  to  the 
profile  at  the  reflecting  region    [j,k],   the  magnetic  azimuth,     cpa  ,   the 
magnetic  dip  angle,     I,  and  the  magnetic  intensity,    \ti\  ,     are  required  to 
calculate  a  reflection  coefficient  matrix    T, 

T  _  r    Tee     Tem  -1  ^  (13) 

T  T 

me  mm 

where, 

T=  T(j,k)     .  (14) 
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Quite  frequently,   a  particular  investigation  of  low  frequency- 
propagation  will  require  only  average  parameters  over  the  entire 
propagation  path.     Thus,   the  magnetic  field  parameters  at  the  center  of 
a  long  propagation  path  in  the  daytime  at  VLF  usually  suffice  to  describe 
the  propagated  field.     An  example  of  a  26  kHz  field  propagated  at  a 
magnetic  azimuth,    cpa  ~  147*     from  Hawaii  will  require  magnetic  para- 
meters given  in  table  1    for  each  member  of  the  set    Lj>k3,     for  a  distance 
of  6180  km  between  transmitter  and  observer. 


Table  1.     Magnetic  Parameters  for  Figures  16  and  17 

[j,k] 


1, 

1 

2, 

1 

2, 

2 

3, 

1 

3, 

2 

3, 

3 

4, 

1 

4, 

2 

4, 

3 

4, 

4 

U|A/m 

cpa   degrees 

I,   degrees 

25.8 

147.8 

-   3.1 

26.0 

147.8 

21.2 

28.4 

144.9 

-26.9 

26.  6 

147.0 

28.2 

25.8 

147.8 

-    3.1 

29.7 

147.  2 

-33.8 

27.0 

146.  5 

31.4 

25.  6 

148.2 

9.5 

2  6.  8 

146.7 

-15.  5 

30.4 

142.  3 

-37.0 

The  amplitude  of  the  field    |Er  |     and  the  fields     JEr>  j  | , 
j  =  1,    2,    3,    *  •  •    ,   were  calculated  as  a  function  of  distance  (fig.    16). 
The  normal  daytime  noon  profile  (fig.    8)  was  employed.     The  magnetic 
parameters,   [1,1],   table  1,   at  the  center  of  the  propagation  path  for 
d  =  3090  km  were  used  for  the  entire  path.     In  figure  17,   the  magnetic 
parameters  for  each  reflecting  region,    [j,k],     given  in  table  1,    were 
inserted  into  the  calculation  so  that  a  separate  reflection  coefficient  for 
each  reflecting  region  [j,k]  was  calculated.     As  the  distance  changes 
it  was  assumed  that  the  reflecting  regions  [j,k]  each  remained  the  same, 
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It  is  noted  that  only  minor  changes  in  the  field  occur  as  a  consequence 
of  a  variable  magnetic  field  at  each  reflecting  region.     It  should  be 
noted  that  this  propagation  crosses  the  magnetic  equator  in  such  a  way 
that  the  dip  angle,   I,   table  1,   changes  from  positive  to  negative  along 
the  propagation  path.     We  therefore  conclude  that  the  value  of  magnetic 
field  at  the  center  of  the  propagation  path  may  be  a  good  assumption  for 
daytime  propagation  over  such  a  path.     Thus,   instead  of  calculating 
10  sets  of  reflection  coefficients,    one  set  for  the  center  of  the  propaga- 
tion path  would  suffice  if  accuracy  of  the  order  of  a  decibel  were 
sufficient. 

The  detailed  nature  of  the  terms  of  the  geometric  series,     Er>  j  , 
j  =  0  to  7    illustrate  typical  behaviour  for  such  terms  as  a  function  of 
distance.     Thus,   at  great  distance,     d/ j,     the  amplitude  decays  almost 

exponentially  (linear  semilog  graph).     At  shorter  distances,     d/j, 

i  l 
small  undulations  exist  in  the  curves.     The  total  field,    ^~      ,     does  not 

L 

$=° 

reflect  much  of  this  detail  but  does  show  undulations  as  a  function  of 
distance  resulting  from  phase  interference  between  terms  of  the 
geometric  series. 

The  exponential  and  two -layer  profiles,   Thrane  (1969)  (fig.    5) 
and  the  Bain  and  May  (1967)  two -layer  profile  (fig.    6)  were  used  to 
calculate  the  amplitude  and  the  phase  correction  as  a  function  of  distance 
at  20  kHz,    see  figures  18,    19,    20,    21,    22,   23,    24,   25,  and  2  6.     The 
phase  correction,     cpc  ,   is  defined 

cpQ   =  -Arg  Er    -  kid    .  (15) 

The  amplitude  is  normalized  to  unity  dipole  current  moment.     Thus, 
the  radiated  power,     Pr  ,   is  given  by: 

Pr  =  0.4(10-ls)  w2  (I0<l)2/z0  (16) 

where     Z0   =  377  ohms. 

-14- 


In  the  presence  of  the  exponential  ionosphere  (fig.    18)  the  inter- 
ference pattern  between  the  ground  wave  and  the  first  ionospheric  wave 
is  clearly  superposed  on  the  inverse  distance  decrement  in  the  amplitude. 
Of  course  at  greater  distances  the  higher  order  ionospheric  waves, 
j  =  2,    3,    *  *  *    ,   become  important  and  also  cause  interference  phenomena. 
The  phase  correction  (fig.    20)  also  shows  this  interference  pattern.     An 
interesting  interference  pattern  between  the  first  ionospheric  wave  and 
the  ground  wave  at  short  distance  and  between  the  various  ionospheric 
waves  at  greater  distance  is  shown  by  the  calculations  in  figures  22  and  23, 
These  graphs  are  based  on  the  two -layer  profile  shown  in  figure  5.     In 
particular,   a  deep  minimum  is  evident  near  600  km.   Calculations  for  an- 
other two -layer  model  given  by  Bain  and  May  (19  67)  (fig.  6)  are  shown  in 
figures  28  and  29.   The  ratio  of  the  first  ionospheric  wave,   j  =  1,   to  the 
ground  wave,   j  =  0,  „ 


20  log  1  -g** 


ztX 


r>  o 


is  illustrated  in  figure  24  for  this  two  layer  model.     Also,   the  corres- 
ponding phase  difference, 


ArgEr>1     -ArgEr> 


r,  o      » 


is  given  in  figure  2  6.     The  corresponding  curves  for  the  exponential 
profile  in  figure  5  are  shown  in  figures    25    and    27.      These  curves 
were  calculated  using  sea  water  ground  constants,     0"  =  5,   £%   -  80. 

However,   the  difference  between  the  land,     0=0.  005,   e2   =  1  5  curves 
and  the  sea  water  curves  is  small.     This  can  be  ascertained  by  com- 
paring figures  18  with  19  and  20  with  21.     A  standing  wave,   as  a  function 
of  distance  at  short  distances,    is  shown  in  figure  24  which  is  apparently 
a  consequence  of  a  wave  trapped  by  the  two-layer  model.     It  is  also  clear 
that  the  phenomenon  completely  disappears  in  the  exponential  model  but 
also  exists  in  the  Thrane  (1969)  two -layer  model  of  figure  5.     It  would  be 
possible  to  observe  such  a  wave  if  an  experimental  method  to  eliminate 
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the  ground  wave  could  be  devised.  Thus,  the  presence  of  the  undulation 
as  a  function  of  distance,  (i.  e.  ,  a  standing  wave),  would  be  an  indication 
of  multiple  D -region  layers. 

The  model  of  the  D -region  illustrated  in  figure  7  is  typical  of  the 
D -region  in  the  United  States.     The  60  kHz  field  as  a  function  of  distance 
was  calculated  for  the  Ft.    Collins,    Colorado,   NBS  transmitter  WWV. 
The  magnetic  azimuth  was  assumed  to  be    cpa   =90°   and  270°  .     The  dip 
angle  was  assumed  to  be    I  =  70°  .     The  magnetic  intensity  was  assumed 
to  be,      \ 34 \   =  40  A/  m.     The  ground  constants    a  =  0.  005  mhos/  m  and 
e2   =  15    were  assumed.     The  two  theoretical  curves,   cpa   =  90°   and  270°, 
are  given  as  solid  lines  in  figures  30,    3l,    32,   and  33.     Measured  data 
of  Morgan  (1966)  are  plotted  for  various  radials  from  the  transmitter. 
The  absolute  field  for  13  kw  radiated  power  is  given  on  the  vertical 
scale  in  decibels  above  one  volt  per  meter.     Figure  30  gives  the  mea- 
sured data  on  a  radial  toward  Palm  Beach,  Florida.     Figure  31  gives 
the  measured  data  on  a  radial  toward  Cape  Fear,  North  Carolina. 
Figure  32  gives  the  measured  data  on  a  radial  toward  Brownsville,   Texas. 
Finally,   figure  33  gives  the  measured  data  on  a  radial  toward  Nantucket, 
Massachusetts.     No  attempt  was  made  to  determine  ground  constants 
more  accurately  or  adjust  profile  so  that  data  would  fit  more  accurately. 
The  agreement  with  theory  is  considered  to  be  quite  excellent.     These 
data  are  quite  remarkable,    since  data  of  amplitude  vs  distance  is  rare 
and  can  be  obtained  only  after  considerable  effort.     One  could  of  course 
modify  the  profile  and  other  parameters  to  obtain  a  more  precise  fit  of 
the  data  for  the  various  propagation  directions,   but  determining  a  para- 
meter or  parameters  to  change  to  obtain  the  desired  result  would 
require  more  research. 

To  gain  some  notion  of  the  extent  of  the  amplitude  perturbation 
that  may  occur  as  the  ground  constants  of  a  propagation  path  are  changed 
or  as  the  ionosphere  changes  with  time,   a  series  of  60  kHz  curves  of 
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amplitude  as  a  function  of  distance  are  presented  in  figures  34  through 
43.     These  curves  are  normalized  to  unity  dipole  current  moment. 
Figures  34,    35,    37,    39,    40,   and  42  employ  the  model  given  in  figure  7. 
Figures  3  6,    38,    41,   and  43  employ  the  model  given  in  figure  8,    i.e.  , 
the  alternate  normal  daytime  noon  model  with  somewhat  less  ionization 
between  40  and  70  km.     This  latter  model  is  considered  to  be  more 
sharply  bounded  than  the  former  (fig.    7).     A  single  nighttime  model 
given  in  figure  9  is  used  on  the  nighttime  curve  on  all  the  figures  34 
through  43. 

Figures  34  and  35  illustrate  propagation  northward,     cpa   =  0° 
and  eastward,     cpa   =  90°   assuming  a  ground  conductivity    0"  =  0.  005 
mhos/m.     Although  both  daytime  and  nighttime  curves  of  amplitude  as 
a  function  of  distance  are  given,   the  principal  change  as  a  result  of  the 
magnetic  azimuth  change  occurs  at  night.     The  change  comprises  a 
shifting  of  the  undulations  such  that  at  a  fixed  distance  both  an  increase 
and  a  decrease  in  amplitude  are  possible  with  a  change  in    cpa  .     The 
numerical  details  of  the  computation  indicate  a  strong  coupling  of  energy 
into  the  whistler  mode.     This  causes  the  phase  of  the  reflection  coeffi- 
cient to  vary  rapidly  with  angle  of  incidence  on  the  ionosphere,    i.  e.  , 
with  distance.     This  phenomenon,    together  with  the  normal  interference 
between  the  ground  wave  and  the  ionospheric  waves  generates  a  compli- 
cated curve  as  a  function  of  distance.     Comparison  of  figure   35  with 
figure  3  6  shows  only  a  change  in  the  daytime  curve.      In  figure  3  6 
we  have  introduced  the  alternate,    more  sharply  bounded  daytime  model, 
figure  8,    with  the  lower  electron  number  density  between  40  and  70  km. 
Although  this  model  does  not  fit  the  data  shown  in  figure  30  as  well,    it 
is,  non  the  less  a  physically  possible  model.    This  gives  us  some  quantitative 
intuition  concerning  the  magnitude  of  perturbations  possible  at  this 
frequency  as  the  lower  D-region  changes  from  high  to  lower  electron 
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density.     Figure  37  should  be  compared  with  figure  34  since  the  only- 
change  between  these  figures  is  an  increase  in  the  ground  conductivity 
from    a  =  0.005    to    o  =  0.  01  mhos/m.     Similarily,    figure  38  should  be 
compared  with  figure  36    (cpa  =  90°,   alternate  daytime  noon  ionosphere 
model  fig.    8).  Figure  28  should  be  compared  with  figure  35.     Again, 

if  the  conductivity  is  lowered  from    0  =  0.  01    to    o  =  0.  001  mhos/m, 
compare  figure  37  with  figure  40  and  figure  3  6  with  figures  43  (o"  =  0.  005  to 
o  =  0.  001).     It  is  apparent  that  one  cannot  guarantee  that  a  particular 
observer  may  not  experience  as  much  as  1  5  dB  decrease  or  increase  in 
signal  at  a  particular  locality  as  a  result  of  ionosphere  temporal  changes. 
Thus,    if  a  disturbance  of  the  ionosphere  occurs  in  the  region     R    in 
figure  12,    the  first  ionospheric  wave,     Er>1,     will  be  affected.     Also 
(figure  14)  the  third  ionospheric  wave,     Erj2   ,     will  also  be  affected.     It 
is  then  probable,    if  the  disturbance  is  local  that  the  second  ionospheric 
wave,     Erj2,  (figure  l3)and  the  fourth  ionospheric  wave,     Erj4,   will  be 
unaffected.     This  disturbance  may  be  a  change  in  the  profile  or  a  lower- 
ing of  the  reflecting  region.     Then,   wave  interference  between  the  dis- 
turbed and  the  undisturbed  ionospheric  waves  could  also  produce  fading 
similar  to  that  described  above  at  60  kHz  in  which  the  entire  propagation 
path  was  disturbed.     The  computer  program  described  in  section  7  will 
accommodate  separate  reflection  coefficients  for  each  of  the  reflecting 
regions,     [j,k].    Also,  a  large  number  of  such  disturbances    Lj»k]     can 
be  accommodated  depending  upon  the  dimension  size  of  the    j,k  matrices. 
Under  such  circumstances,   the  principal  interest  is  amplitude  and  phase 
as  a  function  of  time  rather  than  distance.     Numerical  studies  of  this 
type  are  reserved  for  future  work. 
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3.     THEORY  OF  n-WAVES 

The  method  of  summation  of  the  series  (2)  although  surprisingly- 
simple  tends  to  mask  the  physical  significance  of  the  n-wave  s.     This  is 
a  consequence  of  the  fact  that  the  factor,     G(8,n),    contains  the  Legendre 
function,     Pn(cos8),    which  is  a  real  number  with  alternating  sign.     It 
would  appear  therefore  that  each  n-wave  propagates  only  in  the  radial 
direction.     It  would  also  appear  that  each  n-wave  is  a  standing  wave  in 
the  direction  of  increasing  distance,    d  =  a8.     However  such  a  standing 
wave  must  be  composed  of  traveling  waves  in  the    8 -direction  super- 
imposed upon  each  other,   but  traveling  in  different  if  not  opposite 
directions. 

Using  results  obtained  by  Hob  son  [1896]    and  Barnes  [1907], 

Pn(-cos  9)  =  Pn(cos8)  exp  (T  inn) Qn(cos  8)  sinnn  .  (17) 

This  can  be  written  for  integral    n    only  as 

Pn(cos8)  =  (-l)ftPn(-cos8)  =  (-l)RPn  [cos(tt  -8)].  (18) 

It  is  found  from  this  that 


Pn(cos8)    ~  ("1)n {(cosa  +  isin^)  +  (cosQ  -isinQ)}  (19) 


V  2rr(n+l)  sin  8 


where 


If,* 


TT 


n  =  (n+  J)  (tt-9)  -ii     .  (20) 


d  =  ge  ,  (2i) 


*See  section   6  for  definitions  of  symbols. 
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the  lateral  displacement  of  the  wave  along  the  ionosphere, 

1 


and 


n+f 

s   — 

kog 

K 

c 

(22) 


(23) 


then,   using  (19), 


Pfi(cos0)~  V  i(-l)M-l)^_  (exp(_ikoDs)  +  iexp(ik0Ds)},      (24) 


where 


s]  2TT(n+l)  sin9 


D7  =  D  -  2ng  , 


the  negative  or  reverse  net  lateral  displacement  (distance)  at  the  ionosphere 
around  the  world  the  long  way. 
Hence,    (2)  can  be  written  as 

_  _  M-o  c       Io*>       V       n(n+l)(2n+l)  cfi  Cw    (1  +  RJ3       j  r 

E-3  -"sTT  F^  ^  /7TTf  +n    •  5  Po    ' 

-i         ^=0  V  2n  (n+1)  sin  9 


i    (-1)*+1    [exp( -ik0Ds)  +  iexp(ik0Ds)  }     .  (25) 


The  factor  involving  the  exponential  functions  is  the  sum  of  two  plane 
waves  for  either   outgoing    £ia     or  ingoing    £ia     type  waves  or  rays  tan> 
gent  to  the  harmonic  sphere  of  radius    rn .     Actually  two  sets  of  ortho- 
gonal rays  tangent  to  a  harmonic  sphere,     rfi ,    exist.     Thus,   in  the 
meridian  plane  of  the  terrestrial  sphere, 

B  =  £±1   =  %.     ,  (26) 

or 


n+  2     _ 
k0g 

g 

9 

n+  2 

k0 

• 

ra  =  -^       .  (27) 
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According  to  (10),     s  «  sincpj  .     Thus,    in  general,    there  exist  two 

,(1,2)  f  Pa(-cos9)        iQfl(  -cos  9)^, 
orthogonal  waves    C,  I ±  — w_j.)»     or  four  type 

waves  (corresponding  to  superscripts   l,2,and+,   -signs)  if    Qn     is  the 
second  solution  ( Abramowitz  and  Stegun,    1964)  of  Legendre'  s  differen- 
tial equation.     As    n    increases,     rn     becomes  greater  than    a  +  h,    cpx 
and    sincpi    become  complex,   and    Re  sincpj     becomes  greater  than  unity, 
i.e.,    the  waves  become  inhomogeneous. 

It  is  interesting  to  note  that  (25)  can  be  written  as 

00  00 

Erj  j  =  A   Y    B(n)  exp  (-ik0Ds)  +  iA   V     B(n)  exp(ik0Ds)  .         (28) 

n=°  ft=o 


More  generally  (25)  or  (2)  can  be  written  as 

F  -    C    S       CM    l     1\»    T    P^-C°S9)  *     Qn(-CQS9)    "1 

Er,  I    -  C    I    G(n)  (-1)-    |_  —^ i  -  j 

+  C^G(n)(-l)»    [   P-'-c2°s9'  +iQ,(-coB6)    "j     _      (29) 


a  =o 


The  quantities    A,   B(n),    C,   and    G(n)    are  defined  by  reference  to  (25) 
and  (2).     In    equations     (25)    and    (29),    the  second  terms,     involving 
a  summation  over    n,     represent  the  wave  traveling  around  the  terres- 
trial sphere  in  the  reverse     (-9)    direction.     This  field,    is,    of  course, 

negligible,    except  near  the  antipode  of  the  transmitter  at  low  frequencies. 
Equation  (2)  is  therefore  valid  near  the  antipode  and  should  produce  a 
standing  wave  as  a  consequence  of  the  wave  in  the  reverse  direction. 
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At  short  distances,   as  a  consequence  primarily  of  the  great  lateral  dis- 
placement of  the  n-wave  depicted  in  figure  11  and  the  resultant  severe 
attenuation,    the  second  summation  term  can  usually  be  neglected. 
Multiple  forward    (+8)    and  reverse    (-9)    waves  can  be  readily 
demonstrated  with  the  aid  of  Hobson'  s  ( 1896)    convergent  series  for 

TT  5TT 

T   <  0  <    — r      ;     here, 
6  5 


r(n)  cos{(n+|)9-^  }  i  cosl(n+f)8-lg} 

P»(cose)  =  n    F(ZTi)  L        -    .    "  1  +  2  •  2n+  3  ,.    .    „,3/2 

(2  sinB)  (2  sin  8) 


1  cosl(n+|e    -gj        ^ 

+  2-  4-  2n+ 3-  2n+ 5  ,„     .     ..3/2  +#"J  (30) 

(2  sin  9) 

and 

«,      «■    r-     r(,)r-"nt(n+i)6-i1  !        ^t(n+|)e-^} 

Qn(cosB)=Vn    __L___ — z^teiT         .,    .    „,3/2 

2'  (2sm8)2  (2sm8) 

sin  {(n+ |)8   -^} 

-    -  ••'         ,  (31) 


2-4-  2n+3-  2n+  5  7,     .     ..3/2  J 

(2  sm8) 


which  can  be  combined  as 


,    /  ax  rw  ax  r/  \ '       -    exp  [  Ti(n+i)(TT-9)  ±  i  j  ] 

^(-cosS)       .   QR(-cos8)  _(    .a      r(n)       j  4 

2  X  tt  l"*'     Hn+i)     I 


r(n+^     L  (2TTsin0)^ 

j  exp[  =F  i(n+  |)(TT-9)  ±  i~^] 

+    2-  2n+3  ~     .    173/2 

(2n  sm9) 

expl  T  i(n  +  |)(TT-9)  ±~ 

+  _ - —  + 

2  •  4-  2n  +  3  •  2n  +   5  .„         .     AX5/2 

(2n  sin 9) 
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}.  '32) 


where 

T(n)  1  1 

It  is  apparent  that  the  equation  (32)  splits  the  series  (29)  into  waves 
traveling  in  the  positive  (  +  9)    or  negative     (-0)    direction;    the    +iQfl 
refers  to  waves  in  the     -0     direction  and  the     -  iQn     refers  to  waves 
in  the    +  6     direction.     The  higher  order  terms  of  the  series  (32)  can 
usually  be  neglected  if    n    is  large  or    n~  kia;     the  remaining  terms 
comprise  the  important  terms  for  the  radio  wave  propagation  problem 
under  discussion. 

Equation  (32)  can  be  written, 


Pn(-cos9)  .    Qn(-cos9)    = 

2  a  TT 


("1}     HnTTT  lr,      .  fi-,J   expL  "lk°  Sl  D  J 

v        '  [2n  sin9  J  ^ 


+  i  r  D  "I 

+  exp  I    -ik0   s2        j 

(2.2n+3)[2rrsin9]J/^ 


+   i         2 


5/2 

(2.  4.  2n+3.  2n+5)  [2n  sin9] 


exp 


ik0    s3   £    j    +  •••  }     ,  (33) 
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where  , 

n+z 

-I 

S2     =   — i 

k0g 

S3  =  — i 


Thus,   the  value  given  for     s    approximately  in  (22)  is  a  real  number. 
The  angle  of  incidence  of  a  single  plane  wave  can  be  replaced  more 
precisely  with  a  sum  of  plane  waves  with  angles  of  incidence    si    ,    s2    , 
s3    •  •  •    .     In  short,      s  =  sincpj     is  a  complex  angle  of  incidence.     This  is 
consistent  with  the  general  theory  complex  angle  for  the  n-wave  (10). 

4.     CONCLUSIONS 

Flexibility,    computation  accuracy  and  precision  control  and 
a  quite  general  propagation  theory  for  VLF,    LF  and  MF  terrestrial 
radio  wave  can  be  achieved  with  the  use  of  spherical  wave  functions  of 
integer  order  -when  such  functions  are  used  as  building  blocks  for 
computer  programs.     The  analysis  technique  suggests  a  wide  variety  of 
scientific  applications  to  the  study  of  VLF,    LF  and  MF  ionospheric 
wave  propagation,   only  a  very  few  of  which  have  been  illustrated  in  this 
paper. 

It  has  been  noted  that  wave  interference  and  terrestrial  magnetic 
field  phenomenon  and  profile  electron  density  changes  can  affect  the 
amplitude  of  the  60  kHz  signal  by  as  much  as  15  db  at  a  fixed  distance, 
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especially  at  night.     It  has  also  been  found  that  the  26  kHz  signal  is 
comparatively  insensitive  in  the  daytime  to  the  variations  in  the  magnetic 
dip  angle  when  propagation  occurs  across  the  magnetic  equator.     Thus, 
the  value  of  the  magnetic  field  vector  at  the  center  of  the  path  is  a  good 
approximation. 
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6.     GLOSSARY 

Er   =  vertical  electric  field  in  spherical  system    (r,9,cp)  r  =  radial 

direction. 
Er   0    =  ground  wave  field. 
Er    .    =  ionospheric  wave  field,     j  =  1,    2,    3  •  •  • 

B  =  constant.       B"=  \ig    !£-, 

[±Q   =  permeability  of  space,     p.0   =  4tt(10"7  )  Henry/  m. 
IQi>  -  source  dipole  current  moment,   ampere -meters. 
k  -    =  wave  number  of  air, 

(U  ID 

k  ,    =  -Tix    ~    -    =  k0  (Tix  ~  1.0001  to  1.0003). 
-1        c  c 

c  =  speed  of  light,     c  -  2.  997925(108)m/  s 
f  =  frequency,   Hertz,     f  =  W)/  2tt   . 
a  =  radius  of  the  earth,     a  ~  6.  3  6739(1 0s  )m. 
k_2   =  wave  number  of  the  ground, 


i  m  •     ° 

-2       c   N  e0uo 

e2   =  dielectric  constant  of  the  ground. 

0"  =  ground  conductivity,   mhos/m 

e0   =  permittivity  of  space,     e0  ~  8.  8541 9(10"12) 

G(0,n)  =  n(n+  l)(2n  +   1)    Pn(cos9) 

d  =  distance  along  a  great  circle,     d  =  a6 

(n)  =  the  set  of  integers,     n  =  0,    1,    2,    3  •  •  •    . 

Pn(z)  =  the  Legendre  function,   which  obeys  the  recurrence  relation: 

PM-i(*)=-^71    zPn(z)  -JL-  PB.i(z)       , 
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where 


P,(a)  =  1 
Pi(z)  =  z. 


tn 


k_ 


Re  = 


-la      k_ 


^ln' 


_2a 


Po 


r(l)  r        r(2) 

r(s)    J  L    r(i)  J 


C^a    =   dX)  (k_xa) 


£(-h    =cia)  (k-xa) 


tog        -  ^»      (k0g) 


-2  a 


ilia(k.2a)    ~    IC£>a 


Ci1,s)(z)=7=f     H<V?(z) 


(l    2) 

Ha  '      (z)  =  the  Hankel  function  of  argument  z, 

order    n,    of  the  first  or  second  kind 
(Abramowitz  and  Stegun,    1964). 


(!)#»_,     r(2) 


(z)  =  1 1_  cr;  (z)  +  cj ;  (Z)  j 


(l,2) 

Ca        (z)     satisfies  the  recurrence   relation: 

M1'8)/    ^  2tl+1  (1,2)  (1,8) 

Cft+i    (z) Cn        (z)    -&a_i  (z) 
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where 


if 


d1,2)(z)  -    +     iexp(+iz)    , 


C(!x'2)(z)  =  exp  (+  iz) 


f  =  cf,2)(kr)Pn(cose)         , 


(V2  +  k2 )  f  =  0 

■s     r 

is  satisfied  in  spherical  coordinates,   when    — — =  0. 

/    C3-*2)  (i>2) 

^n    C-ia     =  logarithmic  derivative  of    £-ia 


b-la 


(i,2)'_  a      (1,2)      r  a    r(i,2),  .  i 

c-ia  -  ai  c-ia  -  L~zc*    (z)  J, 


N+    =  positive  ion  number  density. 

N_  =  negative  ion  number  density. 

Ne    =  electron  number  density. 

m+  =  mass  of  the  ion. 

me  =  mass  of  the  electron. 

e      =  electric  charge. 

V      =  collision  frequency,    sec"    . 

In'  \|i-ia    -  t-=^-  In    ^-2a 


R„    = 


--In     C-ia    +  ; ln    Y-s 


Cj   =  effective  reflection  coefficient  defined  by  the 
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matrix  equation: 


where 

t=  r Tee  Te,a"]  , 

T  T 

-"-Bis  mm 


G 


_  r  R9      on  r  1       On 

9  "    L  0      -1  J  '      u*  '"  Po  L  0     -R.  J 


Teo  =  reflection  coefficient  (Johler  and  Harper,  1962) 
for  vertical  electric  polarization  of  the  incident 
and  reflected  ionosphere  wave, 

Tmm    =     reflection  coefficient  (Johler  and  Harper,    1962) 
for  vertical  magnetic  polarization  of  the  incident 
and  reflected  ionosphere  wave, 

T0m     =    conversion  coefficient  (Johler  and  Harper,    1962) 
from  vertical  electric  incident  to  vertical 
magnetic  reflected  ionosphere  wave, 

Tm8     =     conversion  coefficient  (Johler  and  Harper,    1962) 
from  vertical  magnetic  incident  to  vertical 
electric  reflected  ionosphere  wave, 

cpt        =    angle  of  incidence  on  the  ionosphere, 


sincp,  =[-iCfI!   (k'*>l, 

Ci'  (Kg) 


g       '=    radius  to  lower  ionosphere  boundary,   g  =  a  +  h, 

h         =    height  of  lower  ionosphere  boundary  above  the  ground, 

Pr      =    radiated  power, 

Pr   =  0.4(10-12)uo2  (l0l)2  /z0    , 
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cpa     =       magnetic  azimuth,   measured  clock-wise  from 

magnetic  north,    degrees  or  radians, 
|3i|    =     magnitude  of  the  terrestrial  magnetic  field  vector    W  , 
I       =     magnetic  dip  or  inclination,   measured  downward 

from  the  horizontal. 


,     1 

n+  2 

s 

k0g 

r 

n+  2 

1  n 

k0 

k„ 

r>-0 

c 

d      =  ge 

Qn(z)  =   the  second  solution  of  Legendre*  s  differential  equation 
(Abramowitz  and  Stegun,    19  64). 
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Figure  3. 


Amplitude  of  the  ionospheric  waves,   j  =  1,    2,    3  •  •  •    , 
and  the  ground  wave,  j  =  0,   as  a  function  of  distance; 
waveguide  height  =   65  km;    Te0   =   -1,   0"  =  0.  005  mhos/m; 
f  =  10  kHz. 
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Figure  4.  Amplitude  of  the  ionospheric  waves,   j  =  1,    2,    3  ••• 

and  the  ground  wave,  j  =  0,   as  a  function  of  distance; 
waveguide  height  =   65  km;    Tee   =   -1;  a  -  0.  005mhos/  m; 
e2  =15;   f  =  26  kHz. 


-37- 


Oli 


_D   jo   ul6o| 


E 


a> 

X 


o> 


•  H 
CO 

U 


so 

ej 

^  •■-| 

Pi 

o 

a. 

to 

o 

u 


0) 
^     a) 

rH      CXI 

I    o 

I- 


(0 

u 


CO 

rt    ft  d 

v  2 


cj 


«  .2 

O   t3 


C3 

§3 


if} 


a> 

<D 

X 

M 

2 

CXl 

•  H 

h 

.aio/ N    Bo| 


-38 


g 

o 


80 
Height,   km 


Figure   6.         Alternate  two-layer  model  electron  number 

density  profile  of  the  ionosphere  together  with 
the  corresponding    d-c  conductivity,   mhos/m. 
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Figure  7.         Normal  daytime  (noon)  model  electron  number 
density  profile  of  the  ionosphere  together  with 
the  positive  ion  number  density  and  the 
corresponding  d-c  conductivity  resulting  from 
ions  (positive  and  negative)  and  electrons. 
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Figure  8.         Alternate  normal  daytime  (noon)  model  electron 
number  density  profile  of  the  ionosphere  together 
with  the  positive  ion  number  density  and  the 
corresponding  d-c  conductivity  resulting  from 
ions  (positive  and  negative)  and  electrons. 
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Figure  11.        Diagramatic  representation  of  an  n-wave  for  j  =  1  and 

j  =  2,    i.  e.  ,    for  single  and  double  ionospheric  reflections, 
illustrating  principal  regions  of  ground  and  ionosphere 
illumination.     The  geometric -optical  ray  limit  is  dotted. 
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Figure  18.      Theoretical  amplitude  of  the  total  field 
as  a  function  of  distance  using  the 
exponential  model  profile,    figure  5; 

f  =  20  kHz;    o=0.  005  mhos/m;    e2   =15; 

cpa   =  90°;   I  =  70°  ;    |w|   =  36  A/ m. 
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Figure  19.      Theoretical  amplitude  of  the  total  field 
as  a  function  of  distance  using  the  two- 
layer  model  profile,  figure  5;    f  =  20  kHz; 
a=5  mhos/m;    e2   =  80;   cpa   =  90°; 
I  =  70°;    |jt|   =  36A/m. 
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point. 
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Figure  20.       Theoretical  phase  correction  of  the  total 

field  as  a  function  of  distance  using  the 

exponential  model  profile,  figure  5; 

f  =  20  kHz;    o  =  0.  005  mhos/m;    e2   =  15; 

cpa   =  90°;    I  =  70°;    \%\   =  36A/m. 
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Theoretical  phase  correction  of  the  total 
field  as  a  function  of  distance  using  the 
two -layer  model  profile,  figure  5; 
f  =  20  kHz;   a  =  5  mhos/m;    e2   =  80; 
cpa    =90°;   I  =  70°;    \%\   =  36  A/ m. 
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Figure  22.  Theoretical  amplitude  of  the  total  field 
as  a  function  of  distance  using  the  two- 
layer  model  profile,  figure  5;    f  =  20  kHz; 

0  =  0.  005  mhos/  m;    ea    =  1  5;    cpa   =  90°  ; 

1  =  70°;    \U\   =  36A/m. 
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Theoretical  phase  correction  of  the  total 
field  as  a  function  of  distance  using  the 
two-layer  model  profile,   figure  5;  f  =  20  kHz; 

0  =  0.  005  mhos/  m;    e2   =  1  5;    cpa   =  90°  ; 

1  =  70°;    |Ml   =  36  A/m. 
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Figure  24.        Theoretical  amplitude  of  the  first  ionospheric 
wave  relative  to  the  ground  wave  for  the  two- 
layer  model,   figure  5;    f  =  20  kHz; 
a  =  0.  005  mhos/m;   e2   =15.     For  total 
field  see  figure  22. 
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Figure  25.      Theoretical  amplitude  of  the  first  ionospheric 

wave  relative  to  the  ground  wave  for  the 

exponential  model,   figure  5;   f  =  20  kHz; 

a  =  5  mhos/m;   e2   =  80.     For  total  field 

see  figure  18. 
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Figure  2  6.       Theoretical  phase  of  the  first  ionospheric 
wave  relative  to  the  ground  wave  for  the 
two -layer  model,   figure  5;    f  =  20  kHz; 
o  =  0.  005  mhos/ m;    e2   =  15.     For  total 
field  see  figure  23. 
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Figure  27.       Theoretical  phase  of  the  first  ionospheric 
wave  relative  to  the  ground  wave  for  the 
exponential  model,   figure  5;    f  =  20  kHz; 
0=5  mhos/m;    e2   =15.     For  total  field 
see  figure  20.        ,  _ 
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Figure  28. 
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Theoretical  amplitude  of  the  local  field  as 
a  function  of  distance  using  the  alternate 
two -layer  model  profile,    figure   6;    f  =  20  kHz; 
a  =  0.  005  mhos/m;    e2   =15;   cpa   =  90a  ;    1=  70°  ; 
|M|   =  36A/m. 
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Figure  29.       Theoretical  phase  correction  of  the  total 
field  as  a  function  of  distance  using  the 
alternate  two -layer  model  profile,  figure  6; 
f  =  20  kHz;   a  =  0.005  mhos/m;    e2   =  15; 
cpa   =  90°;   I  =  70°;    \u\   =  36A/m. 


-62- 


1 00|?- 


ax 

?    50< 
a. 


60- 


ao«- 


30- 


20«- 


10- 


Theoretica 
<*>a=270° 


0.0   500.0   1000.0 


1500.0  2000.0 
D1ST.  KM 


2500.0      3000.0 


Figure  30.        Theoretical  amplitude  of  the  total  field  as  a 
function  of  distance  for  the  NBS  Ft.    Collins, 
Colorado  transmitter  (normalized  to  13  kw 
radiated  power)  using  the  normal  daytime 
noon  profile  model,   figure  7;    £  =   60  kHz 
a  =  0.  005  mhos/m;    e2   =  1 5;    cpa   =  90°  ;    and 
270°;   1=70°;    |tt|=40A/m.   Data  (Morgan, 
1966)  measured  on  a  radial  from  the 
transmitter  toward  Palm  Beach,   Florida 
are  also  shown. 
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Figure  31.        Theoretical  amplitude  of  the  total  field  as  a 
function  of  distance  for  the  NBS  Ft.    Collins, 
Colorado  transmitter  (normalized  to  13  kw 
radiated  power)  using  the  normal  daytime 
noon  profile  model,   figure  7;    f  =   60  kHz; 
a  =  0.  005  mhos/m;    e2   =  1 5;   cpa  =  90u  and  270* 
1=  70°;    |tf|   =  40  A /m.     Data  (Morgan, 
1966)  measured  on  a  radial  from  the 
transmitter  toward  Cape  Fear,   N.    C. 
are  also  shown. 
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Theoretical  amplitude  of  the  total  field  as  a 
function  of  distance  for  the  NBS  Ft.    Collins, 
Colorado  transmitter  (normalized  to  13  kw 
radiated  power)  using  the  normal  daytime 
noon  profile  model,  figure  7;    f  =   60  kHz; 
a  =  0.005  mhos/m;    e2   =  15;   cpa   =  90°  and  270s 
I  =  70°;    \U\   =  40  A /m.     Data  (Morgan, 
1966)  measured  on  a  radial  from  the 
transmitter  toward  Brownsville,    Texas 
are  also  shown. 


-65- 


100 


60- 


en 

?    50 
a. 

z: 


40- 


30- 


20'- 


10- 


Datum   Point 


0.0        500.0      1000.0 


1500.0     2000.0 
D1ST.  KM 


2500.0      3000.0 


Figure  33.       Theoretical  amplitude  of  the  total  field  as  a 
function  of  distance  for  the  NBS  Ft.    Collins, 
Colorado  transmitter  (normalized  to  13  kw 
radiated  power)  using  the  normal  daytime 
noon  profile  model,   figure  7;    f  =  60  kHz; 
a  =  0.005  mhos/m;   e2  =  15;   cpa  =  90°  and  270" 
I  =  70°;    |tt|  =  40  A/ m.     Data  (Morgan, 
1966)  -measured  on  a  radial  from  the 
transmitter  toward  Nantucket,   R.   I. 
are  also  shown. 
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Figure  34.        Theoretical  amplitude  of  the  total  field  as  a 

function  of  distance  using  the  normal  daytime 
noon  profile  model  figure  7;    and  night- 
time model,   figure  9;    f  =  60  kHz; 
a  =  0.  005  mhos/m;    e2   =15;   cpa   =  0*  ; 
I  =  70°;    \U\   =  40  A/m. 
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Figure  3  5.        Theoretical  amplitude  of  the  total  field  as  a 

function  of  distance  using  the  normal  daytime 
noon  profile  model  figure  7;   and  the  nighttime 
model,  figure  9;    f  =  60  kHz;    a  =  0.005  mhos/m; 
e2   =  15;   cpa   =  90°;   I  =  70°;    \u\   =  40  A/ m. 
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Figure  36.        Theoretical  amplitude  of  the  total  field,   as  a 

function  of  distance  using  the  alternate  normal 
daytime  noon  profile  model,  figures  8,    9  and 
the  nighttime  model,   figure*?;     f  =  60  kHz; 
O  =  0.  005  mhos/m;    e2   =  15;   cpa   =  90°;    I  =  70°  ; 
\U\   =  40  A/m. 
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Figure  37.       Theoretical  amplitude  of  the  total  field,  as  a 
function  of  distance  using  the  normal  daytime 
noon  profile  model,    figure  7;   and  the  nighttime 
model,  figure  9;    f  =  60  kHz;   a  =  0.  01  mhos/m; 
e2   =  15;    cpa   =  0°;   I  =  70°;    \u\   =  40  A/ m. 
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Figure  38.       Theoretical  amplitude  of  the  total  field  as  a 

function  of  distance  using  the  alternate  normal 
daytime  noon  profile  model,  figures  8,    9  and 
the  nighttime  model,   figure  9;    f  =   60  kHz; 
a  =  0.  01  mhos  /  m;    e2    =  1  5;    cpa   =  0°  ; 
I  =  70°;    |w|   =  40  A/m. 
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Figure  39.       Theoretical  amplitude  of  the  total  field  as  a 

function  of  distance  using  the  normal  daytime 
noon  profile  model,  figure  7  and  the  nighttime 
model,  figure  9;    f  =  60  kHz;    0"  =  0.  01  mhos/m; 
e2   =  15;    cpa   =  90°;    I  =  70°;    |m|   =  40  A/m. 
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Figure  40.       Theoretical  amplitude  of  the  total  field  as  a 

function  of  distance  using  the  normal  daytime 
noon  profile  model,   figure  7  and  the  nighttime 
model,  figure  9;    f  =   60  kHz;    a  =  0.001  mhos/m; 
e2    =15;    cpa   =  0°  ;   1  =  70°;    \  W  |   =  40  A  /  m. 
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Figure  41.        Theoretical  amplitude  of  the  total  field  as  a 

function  of  distance  using  the  alternate  normal 
daytime  noon  profile  model,   figures  8,    9  and 
the  nighttime  model,  figure  9;    f  =  60  kHz; 
a  =  0.  001  mhos/m;    e2   =15;   cpt    =  0°  ; 
I  =  70°;    \u\    =  40  A/m. 
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Figure  42.        Theoretical  amplitude  of  the  total  field  as  a 

function  of  distance  using  the  normal  daytime 
noon  profile  model,    figure  7  and  the  nighttime 
model,  figure  9;    f  =  60  kHz;    o  =  0.001  mhos/m; 
e2   =  15;    cpa   =  90°;    I  =  70°;    |g|   =  40  A/m. 
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Figure  43.        Theoretical  amplitude  of  the  total  field  as  a 

function  of  distance  using  the  alternate  normal 
daytime  noon  profile  model,  figures  8,    9  and  the 
nighttime  model,   figure  9;    f  =  60  kHz;    a  =0.0 
a  =  0.001  mhos/m;    e2   =  1 5;   cpa   =  90°; 
I  =  70°;    \U\   =  40  A/m. 
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7.     COMPUTER  PROGRAMS 

Computer  programs  RAYGUIDE  and  RAGIDOO  are  shown  in  the 
computation  plan,   figure  44.     A  program  DIONS  necessary  for  the  cal- 
culation of  reflection  coefficients  from  profiles  and  SORTDION,    which 
smooths  the  phase  of  the  reflection  coefficients  and  references  the  phase 
to  the  reflection  height  are  also  presented.     Thus,   DIONS  generates 
data  cards  for  SORTDION,  and  SORTDION  generates  data  cards  for 
RAYGUIDE  or  RAGIDOO.     RAYGUIDE  calculates  four  ionospheric  waves 
or  hops  simultaneously  and  generates  data  cards  for  RAGIDOO,   which 
can  add,    if  necessary,   a  large  number  of  additional  hops  or  ionospheric 
waves.     RAGIDOO  can  start  with  the  zero  order,   j  =  0,   and  calculate 
hops  sequentially.     But  it  is  not  as  fast  as  RAYGUIDE,   which  calculates 
four  hops  simultaneously.     For  the  majority  of  applications  the  first 
four  hops  suffice. 

The  computer  program  SORTDION  was  written  by  Wieder*  (19  67) 
and  the  program  DIONS  was  written  by  Berry**  (19  65).     The  latter  was 
discussed  in  some  detail  recently  by  Johler  (1967).     The  program 
RAYGUIDE  was  written  by  the  authors.     RAGIDOO  is  a  modification  of 
RAYGUIDE  to  calculate  a  large  number  of  terms  of  the  j -series,   which 
was  originally  written  by  Wilson***  (19  67)  and  modified  by  the  authors. 

The  Legendre  function,     Pn(z)     is  calculated  exactly  in  RAYGUIDE 
for  integral    n,   from  the  recursion  formula, 


Pn(z)  =  Un;1)z    PB-i(z)   -^  Pn-2(z)         ,  (33) 

PoU)  =  1  (34) 


n  n 

where 


*      B.    Wieder  (1967),   Private  Communication,   ITS,   ESSA,   Boulder,    Colo. 
**   L.  A.   Berry  (1965),   Private  Communication,   ITS,    ESSA,    Boulder,    Colo. 
***J.   H.    Wilson  (1967),   Private  Communication,   ITS,    ESSA,    Boulder,    Colo. 
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and 

Pi(z)  =  z  .  (35) 


Both  integer  and  non-integer  values  of  n  can  be  used  in  a 
formula  based  on  the  saddle  point  approximation  of  the  integral  of 
LaPlace  (letting    z  =  cos  8)    as 


2  cos  f(n+J)  6   -in 
Pn(cos8)    ~    —  .  (36) 

J  2n  (n+  1)  sin9 


This  formula  becomes  inaccurate  for  sin  8  ~  0  (short  distances, 
d  =  a  8)  and  n  small.  For  n  large  the  alternative  (3  6)  is  available  in 
the  computer  program. 

The  spherical  wave  functions  are  also  calculated  from  the 
recursion  formulas, 

r(i,3)  .  2n  +  1       (12)..       r(i,2)  .    . 

U  +  i    (z)  =  — - —    £n         (z)   -CB-i    (z)     ,  (37) 

z 


starting  with  the  values, 


and 


d1,3)(z)  =  *iexp[  ±iz]  (38) 


G(_V2)U)  =  exp  [  ±  iz  ]     .  (39) 


Using  a  computation  technique  of  Goldstein  and  Thayler  [1959], 
the  computer  calculates  the  spherical  wave  functions  (37),   (38),  and  (39), 
using  subroutines   RHOSBN    and    ZETAIM.     The  ratio  of  the  derivative 
of  a  spherical  wave  function  to  the  function  can  also  be  determined  by- 
recursion  formulas  ( Johler  and  Berry,    1962),    which  yield, 
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c?'">'w  . 

1 

n 

da'2)(z) 

n 
z 

Cn-l      (z) 
Gr-1     (z) 

z 

(40) 


where 

r(i,2)y,    v 

r(l,2)     /      X 

Co  (z) 

For  large  complex  argument,  z  =  k2a,     the  Debye  approximation  can  be 
employed.     This  yields  ^ 


dl}    (he  a)    _      /  n(n  +  1) 
dl}  (kea)  ,/     (k2a)2 


-  1  •  (41) 


(Im  k2  a  <  0) 

This  very  simple  approximation  is  quite  excellent  for    n(n  +  1)    ^  (k2a)2 
But  the  series  of  zonal  harmonics  converges  to  great  accuracy  before 
n(n  +  1)    grows  as  large  as     |(k2a)2|    .     In  fact,   the  series  converges 

rapidly  just  beyond    n~kia.      Alternately,    the  function       n    .    2   .       can 

(1)  (1)'        *R   (ksa) 

be  readily  used  in  (40)  by  replacing  the    Cn       and    Cfi         function  with 

\|iR     and    tyR      functions.     Thus, 


C,      (k2a)  C  (k2a) 

dl}  (k2a)  ^  (k2a) 


(42) 


/TTz 

where     ^)n(z)  =   /  —    JR+ l(z)     and    JB+ jjz)     is  Bessel1  s  function  of  order 

v        Z  2  2 

n  +   \    and  argument,      z  (Abramowitz  and  Stegun,    1964).      Also, 
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—     ill    (z)  =  \li«'(z).     These  functions  (40,    41  and  42)  are  either  calculated 
dz     Tft\  fa 

directly  in  RAYGUIDE  or   function  PSIR  is  used.     The  latter  function 

uses  real  argument  only. 
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Figure  44.       Conceptual  plan  of  Program  RAYGUIDE. 
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